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Abstract 


This report documents the techniques used to filter quantities on a stretched grid gen- 
eral circulation model. Standard high-latitude filtering techniques (e.g. using an FFT 
to decompose and filter unstable harmonics at selected latitudes) applied on a stretched 
grid are shown to produce significant distortions of the prognostic state when used 
to control instabilities near the pole. A new filtering technique is developed which 
accurately accounts for the non-uniform grid by computing the eigenvectors and eigen- 
frequencies associated with the stretching. A filter function, constructed to selectively 
damp those modes whose associated eigenfrequencies exceed some critical value, is used 
to construct a set of grid-spaced weights which are shown to effectively filter without 
distortion. Both offline and GCM experiments are shown using the new filtering tech- 
nique. Finally, a brief examination is also made on the impact of applying the Shapiro 
filter on the stretched grid. 
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1 Introduction 


Version 2 of the Aries/GEOS dynamical core (Suarez and Takacs,1995) was developed at 
the Goddard Space Flight Center for use in data assimilation and climate prediction appli- 
cations. In this paper we present a generalization of this dynamical core that allows for a 
non-uniform latitude/longitude (A <f>j, AA t *) grid. This capability allows use of the model for 
regional data assimilation with the Goddard Earth Observing System (GEOS) Data Assim- 
ilation System (DAS) (see Rood, 1996) and for regional prediction studies using the NASA 
Seasonal to Interannual Prediction Project’s (NSIPP) coupled prediction system. The use 
of this dynamical core in conjunction with a non-uniform grid has been demonstrated by 
Fox-Rabinovitz et at. (1997) for the Held-Suarez (1994) dynamics benchmark. That work 
demonstrated the usefulness of achieving high resolution results over selected regions with 
much greater efficiency than employing high-resolution uniformly over the sphere. 

The Aries/GEOS dynamical core relies on two filtering strategies to ensure stability and 
efficiency. The first, high-latitude filtering, is done to control linear instability due to the 
converging meridions near the pole. This filter selectively damps the tendencies of fastest 
modes. Thus, on a uniform grid, the filter acts near the poles on the smallest scales. The 
second technique is the application of the Shapiro (1970) filter. This filter is applied globally 
to damp small-scale dispersive weaves and to prevent computational nonlinear instability 
(Phillips, 1959) to occur. This report will examine the properties of these two filtering 
techniques on the stretched grid, and develop a stretched grid convolution filter applicable 
to controlling linear instability in high latitudes. 


2 Stretched Grid GEOS GCM 


Following the work of Fox-Rabinovitz et al. (1997), the stretched grid dynamical core 
was interfaced with the uniform grid physics package of the GEOS General Circulation 
Model (GCM). While the ultimate goal will be to perform the physics on the stretched grid 
itself, this more conservative approach was used to avoid unforseen problems within the 
parameterizations in regions of high resolution. In practice, a copy of the dynamic state 
variables (u, u, 6 , q , p s ) are transformed (via cubic interpolation) to the uniform grid to force 
the physics packages. The tendencies of the dynamic state variables due to the physics are 
then transformed to the dynamics stretched grid to act as external diabatic forcing in the 
time integration scheme. As a result, the dynamics state is always updated and preserved 
on the stretched grid. 

The generation of the stretched grid used for this study employs a cosine mapping function 
to smoothly create grid increments which vary with wavenumber 1 in the longitudinal 
direction, and increase to a maximum at the farthest pole in the latitudinal direction. 
These longitude and latitude grid increments are shown in Fig. 1. For this study, the 
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uniform grid resolution is defined by a ‘2.5°x2.0° lon-lat. grid (144x91 points). These points 
are re-distributed in longitude and latitude so as to produce uniform 1.0°xl.0° resolution 
over the eastern United States (~ 28N-48N, 100W-70W). Figure 2 shows the resulting 2- 
dimensional grid structure over the globe. Note that in the coarsest region (the South-East 
quadrant) an effective resolution of ~ 5°x4° is used. 

Starting from an arbitrary, uniform grid January simulation, initial conditions were trans- 
formed and balanced on the stretched grid to begin stretched grid forecasts. The top panel 
of Fig. 3 shows the initial zonal mean temperature distribution and the day 5 temperature 
change from the intitial condition. During this 5-day forecast, a stratospheric warming 
event has taken place (between day 4 and 5) which is associated with strong cross-polar 
flow. The bottom panel shows the temperature change and day 5 zonal wind at 0.2 mb 
from latitudes 50N to 90N (longitude 0 at bottom). We see that near the poles in the region 
of highest longitudinal resolution, computational small scale features have developed. This 
occurs even though the timestep chosen for the integration was adequately reduced to reflect 
the increased resolution in the area of interest. It was noted in Fox-Rabinovitz et al. (1997) 
that, even for the simple- physics Held-Suarez test, excessive polar noise developed when 
total stretching factors ( ^ m ax ) became large (> 16). In those cases, ad hoc strengthening 
of the polar filters was required. It is clear from these results that, while the standard polar 
filtering technique is effective in preventing linear instability in many cases, it is not robust 
enough for use on the stretched grid for general applications. 


3 Uniform Grid High-Latitude Filter 


To understand why the standard high-latitude filtering technique is not adequate when using 
the stretched grid, we must first review the technique employed on the uniform grid. As 
pointed out in Suarez and Takacs (1995), polar Fourier filters in the Aries/GEOS dynamical 
core are applied to the tendencies of all prognostic variables. This is done to avoid linear 
computational instability due to the convergence of the meridians near the poles. The 
filter acts poleward of a critical latitude 4> c (nominally 45°), and its strength is gradually 
increased toward the pole by increasing the number of affected zonal wavenumbers and the 
amount by which they are damped. 

Consider the linearized, one-dimensional shallow water 
by: 

du g dh 

dt a cos <f> d\ 

dh _ H du 
dt a cos 4> d A 


equations with no mean flow given 

( 1 ) 

( 2 ) 


2 



( 3 ) 


Discretizing in space on a staggered C-grid, we may write 

_ g f hj+i ~hj \ 

dt a cos<f> V AA J 

d_K = H / «,•+! ~ «<-i 
dt a cos<fi ^ AA 

Here, AA = jj= and IM is the zonal dimension. Assuming wave solutions of the form: 

■U l+ 1 (f) = fie‘W + a) AA - ,/ *) , 

fti(t) = he , 



(4) 


(5) 

( 6 ) 


where i denotes the grid-space location in longitude and i = \f^ T, it can be shown that the 
frequency associated with each harmonic component k is given by 


v = ± 


2y5F 

a cos <j> A A 


sin (^^r) • 


(7) 


We see that the frequency increases as a function of wavenumber and as we approach the 
pole. 


Using the leapfrog time scheme, linear stability requires that 

uAt< 1. (8) 


This forces restrictions on the timestep governed by 

x ^ 1 acos(f)AX 1 
-2 TgJr sin(/t^) • 


(9) 


To eliminate the requirement that the timestep goes to zero as we approach the pole, a 
wavenumber-dependent filter, F*, is applied to the time tendencies of the prognostic fields 
resulting in a modified frequency function given by 


v = ± Fk 


2y fffi 
a cos 4> A A 




( 10 ) 


By requiring that the frequency be no larger than i/ MAX (e.g., the unfiltered value associated 
with the shortest wavelength at the critical latitude <f> c ), the functional form of F k may be 
obtained: 


F k = min 1, 


( cos<^> 
cos <f) c 


1 


sin (fc^r) 




an 
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Here, n is an arbitrary factor used by Suarez and Takacs (1995) to increase the strength of 

the filter. While not needed for these linear studies, it has proven beneficial to use n = 2 for 

the non-linear uniform grid GCM. Unless otherwise stated, results shown here were made I 

using n = 1. It is important to remember that applying the filter in this manner simply 

slows down each harmonic component by an amount which produces linear stability. More ? 

precisely, the computational propogation speed of each zonal harmonic is adjusted to keep # 

its transport to no more than 1 AA per timestep. The filter has no impact on the magnitude j 

nor energy of the harmonic. j 

Figures 4 and 5 show the filter function at latitudes 85N and 60N, respectively. For these 

results, the zonal dimension wets set to IM— 144 (corresponding to a 2.5° resolution). Near : 

the pole (85N) the filter strength is quite strong with only wavenumbers 0-5 remaining j 

untouched. However, away from the pole (60N) the filter is much weaker with only moderate j 

damping even at the smallest scales. Also shown in Figs. 4 and 5 are the equivalent j 

grid-space weighting coefficients at various longitudes obtained through convolution of the j 

spectral filter functions. Since the grid is uniform, the weighting stencil is identical for each | 

longitude location. For the strong spectral filter function near the pole, the convolution j 

weights are quite broad in longitude and have a peak amplitude which is relatively small = 

(0.3) . Away from the pole where the filter strength is weak, the convolution weights resemble 

that of a Delta function whose amplitude is near 1.0 at the grid-point location and almost ! 

0. 0 elsewhere. The breadth of the convolution weights are thus intimately connected to the = 

physical scale of the waves being filtered. 1 

To illustrate the effect of filtering, a set of initial conditions consisting of the single harmonics 

1, 3, 5, and 7 (see Fig. 6) is used. To benchmark this test, the integration is first run 
without filters by using a sufficiently small timestep which is linearly stable, Fig. 7. These 
single harmonics simply oscillate in time with no translational motion. The test is then 
repeated using a timestep comparable to that used for GCM simulations. Without filtering, 
instability occurs within a few iterations (not shown). Figure 8 shows the results from the 
uniform grid using the standard FFT polar filter described above. The filtering associated 
with latitude 85N was used for these runs. With tendency filtering, the solutions look quite 
similar to the non-filtered case with no discernable instabilities. Figure 9 compares the 
filtered and non-filtered solutions at longitude 0°, with additional wavenumbers 9 and 12 
also shown. As the wavenumber increases, the frequency is increasingly slowed compared 
with the non-filtered run. Also, the frequencies for wavenumbers 7, 9, and 12 are identical 
for the filtered run since they are forced to be no greater than the maximum allowed value. 

However, the magnitude of the oscillation remains untouched, with no distortion of the 
waveform . 
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The experiment is then repeated using the stretched grid. The uniform 2.5° grid is stretched 
to 1° between longitudes 75E and 105E, Fig. 10, with initial conditions shown in Fig. 11. 
Figure 12 shows the results for the stretched grid using a sufficiently small timestep to 
ensure stability without filtering. The stretched grid solution is very similar to the uniform- 
grid solution, creating negligable distortion due to the stretching. The experiment is again 
repeated using a larger timestep comparable to that used for GCM simulations and thus 
requiring tendency filtering, Fig. 13. It is clear that the use of the standard FFT filter 
has created significant distortions of the pure waveforms and has affected the solution at 
all longitudes. In addition, significant small-scale features have been created even though 
linear stability has been achieved. Arbitrarily increasing the strength of the filter to n = 2 
creates even further distortions, Fig. 14. A snapshot of the solution at 3 hours is shown in 
Fig. 15 clearly demonstrating the inability of the standard FFT filter to adequately control 
noise and prevent distortions on the stretched grid. 


4 Stretched Grid High-Latitude Filter 


Since the standard filtering technique is done in spectral space, it is clear that the functions 
upon which stability is based (the trigonometric zonal harmonics) must have the same 
relevance on the stretched grid for the filtering technique to be valid. Figure 16 depicts the 
single harmonic wavenumber 3 on the uniform grid (top panel) and on the stretched grid 
(center panel). Here, the term “wavenumber 3” is used to describe a field whose structure 
is repeated three times in the longitudinal direction. Using this definition we see that for 
each zonal harmonic k there is an associated geophysical scale whose wavelength L k is given 
by: 


2 7T CZ COS(f> 


This physically-based wavelength is independent of any finite-resolution grid used to con- 
struct the field. Assuming this physically-based wave propogates uniformly in the zonal 
direction, the filter strength required to slow it down must clearly be a function of the 
local zonal grid increment just as it is a function of latitude (i.e., the required strength to 
restrict transport to no more than 1 AA per timestep would need to be greater in an area 
of high resolution than in an area of coarse resolution). However, as shown above for the 
standard filtering technique, there is no longitudinal dependence on the filter strength F k . 
Moreover, since the FFT implicitly assumes a uniform grid distribution, the FFT interprets 
this wavenumber 3 example as a field which is made up of both longer and shorter wave 
components (Fig. 16, bottom panel). The varying strengths of the filter as defined by (11) 
for these wavenumbers would be misapplied if the intent was to filter the geophysical scale 
associated with wavenumber 3. 
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The top panel of Figure 17 shows the FFT harmonic decomposition of the wavenumber 3 
example for both the uniform grid and the stretched grid. Interestingly, the wavenumber 3 
harmonic is almost completely missing from the stretched grid interpretation. Instead we 
see a spread of harmonics to either side of the actual input wavenumber. The bottom panel 
of Fig. 17 repeats this analysis for each wavenumber from 0 to IM/2. Each single harmonic 
wave component is projected onto the strethed grid, and then used as input to the FFT 
analysis. As before, each component is interpreted as having both longer waves and shorter 
waves, with little energy in the actual wavenumber used as input. 

Figure 18 examines the convolution weights obtained from the standard FFT filter but 
applied on the stretched grid. The top panel, as before, shows the spectral filter function 
independent of longitude. The bottom panel of Fig. 18 shows the convolution weights 
plotted on the stretched grid. We see that the peak amplitude of the weights is constant 
in longitude. However, in the fine-resolution region the convolution weights have narrowed 
in longitude while in the coarse-resolution region they have broadened. This implies that 
the effect of the standard filtering technique on the stretched grid will be to filter less in 
the fine-resolution region and filter more in the coarse-resolution region, exactly opposite to 
that required from previous considerations. 

It has become clear that to accurately filter on a stretched grid, we must first determine 
the stretched grid basis functions upon which the model solutions may be projected. We 
again start by writing the discretized linear shallow water equations for the stretched grid: 


du i+± 

9 

dt 

a cos 0 

dhi 

H 

dt 

a cos (j> 


h i + 1 hj 


/ 

H /«,■+! 


Here, is the distance in longitude between mass points, and 


AA - , i + AA 
AA t = ^ '-± 


(13) 

(14) 

(15) 


Using (13) and (14), the wave equation is formed given by: 


d 2 hi gH 


dt 2 [a cos cj >) 2 

Assuming a solution of the form: 


hi hi hi — i 


AA;AA, , i AA;AA,_i 

‘“T 0 1 -> 


2 J 


hi(t) = h,e~^ , 


(16) 


(17) 
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we find 


— i/ 2 (ac os 4>A\)* 

gH 


h{ — 


AA 2 


A \{ A A ’_i_x 


- 2 


+ 


A A 2 


A A • . i_ AA,_ x 

1 ' 2 

AA 2 

AA t -AA. j 


f H + 1 

j h t 

^i-i i 


( 18 ) 


where AA = Note that for the uniform grid (AA, = AA) the weights used in (18) simply 
become 1, -2, 1. Defining: 


AA 2 


r i = 




AA - , i AA- i ’ 
,+ 2 * 2 

AA 2 

AA.AA - , i ’ 


»+5 


AA 2 


AA.AA- i 


(19) 

(20) 
( 21 ) 


we may write (18) in matrix form as: 
/ ~2ri rf 


r 2 — 2r 2 r 2 


’• — 2r IM .j 


h = A , 


( 22 ) 


r IM 


— 2r IM j 


where A is the eigenvalue associated with the eigenfrequency v. 


A = - 


v 2 (a cos 4 > A A) 2 


(23) 


From (22) and (23), w r e may compute the eigenvalues (or eigenfrequencies) and correspond- 
ing eigenvectors associated with the stretched grid (see Appendix for detail). Figure 19 
shows the eigenvalues computed by solving (22) for both the uniform grid and the stretched 
grid, in addition to the analytic uniform grid solution associated with (7). Here we have 
ordered the eigenvalues by magnitude. For the uniform grid, the analytic and numerical 
solutions are indistinguishable. For the stretched grid, we see a significant increase in the 
magnitude of the eigenvalues associated with the highest mode index (smallest scales). 
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The eigenvectors associated with a sampling of modes are shown in Fig. 20 for the uniform 
grid and Fig. 21 for the stretched grid. The uniform grid solutions are simply discretized 
sines and cosines, analogous to the analytic solution. The eigenvectors for the stretched grid 
are similar for the low-index (slow) modes. However, as the mode index (or eigenfrequency) 
increases the structure of the eigenvector is confined to smaller and smaller regions. For the 
highest eigenfrequency, the eigenvector is totally concentrated in the finest resolution area. 


By ordering the eigenvalues and associated eigenvectors in this manner, a filter may be 
constructed to selectively damp those modes whose associated eigenfrequencies are faster 
than some critical value. Using (23) we see that 


, IAI^vW 

V = ± — -r . 

a cos <p A A 


(24) 


As in the uniform grid case, the maximum frequency allowed will be defined as the unfiltered 
frequency associated with the shortest wave at the critical latitude <f> c : 


2^H 

VuAX — / A \ 

Qj cos <p c A^ m in 


(25) 


Requiring that the filtered frequency be no faster than the maximum frequency allowed, 

F v < i'max i ( 26 ) 

the filter function for the stretched grid becomes: 

2 AA cos <f> 


F = min 


1 , 


|A|iAA min cos^ c j 


(27) 


4.1 Convolution Filter 

To construct the stretched grid convolution filter, consider for example the height field 
defined by a vector array h constructed from the set of modes found in Section 4: 

h = Md , (28) 

where M are the stretched grid modes or eigenvectors and A are the mode projection 
amplitudes. The amplitudes may be explicitly obtained by 

A = M _l h . (29) 

A new set of amplitudes are now created by multiplying (29) by a diagonal matrix F defined 
by the convolution filter function (27): 

A f eF 4 = F M -1 h . (30) 
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( 31 ) 


The reconstructed filtered field, h F , is therefore given by: 

h F = M A F = M F M" 1 h . 

Weights in Grid Space 

The weights obtained in grid space are analogous to the weights obtained through convolu- 
tion of the standard FFT for the uniform grid. Through convolution we may now examine 
in spectral space the response of these weights as a function of longitude. The top panel 
of Fig. 22 shows the spectral response of the stretched grid convolution filter for various 
longitude locations, while the bottom panel shows the grid space weighting coefficients. We 
see that the amplitude of the coefficients is now a function of longitude while the breadth of 
the stencil is constant. This is exactly opposite to that obtained by the use of the standard 
FFT filter applied to the stretched grid, with the impact being that the strongest filtering 
now takes place in the region of highest resolution. The results from using the stretched grid 
convolution filter are shown in Fig. 23. The stretched grid convolution filter successfully 
prevents linear instability with no distortion of the waveforms nor generation of additional 
small scale features. 


4.2 GEOS GCM Results 

Having seen the effectiveness of using the stretched grid convolution filter in simple off-line, 
gravity wave test cases, we now proceed to analyze its performance within the GEOS GCM. 
For this test, the forecast experiment described in Section 2 is repeated by replacing the 
standard FFT polar filtering technique with that of the stretched grid convolution filter. 
Figure 24 once again shows the initial zonal mean temperature distribution and the day 
5 temperature change, as well as the 0.2 mb day 5 temperature change and zonal wind 
fields. Compared with Fig. 3, the stratospheric warming event is again simulated but with 
no evidence of small-scale noise near the polar, high-resolution region. It should be noted 
that the noise generated with the standard filtering technique was primarily evident at 
high stratospheric altitudes where the wind field is particular!) 7 strong. At lower levels the 
standard filter and the stretched grid convolution filter produce very similar results. Figure 
25 shows the sea-level pressure and 300-mb height fields at day 5 for the two simulations. 
The two solutions are virtually identical with no evidence of noise in either case. 


4.3 Uniform Grid Interpretation 

While the preceeding sections derived exactly the proper filter function for the linearized 
shallow' water equations on a stretched grid, it is possible to approximate this solution 
through a uniform grid interpretation. Recall from Section 3 that the filter function derived 
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for the uniform grid for n = 1 was given by 


Fk — min 1, 


j COS (j) 

1 )1 

ycos (p c 

sin(*^)/_ 


( 32 ) 


for k = 0, This may be rewritten as: 


for 


F it k' = min 


a cos <f>AAi 

1 )1 

, (33) 

acos^> c AA m j n 

sin 

II 

o 

*13 


(34) 

N- - 2?r 

' A A; • 


(35) 


We see that for the uniform grid, AA; = (32) and (33) are identical. However, for the 

stretched grid we define N s at each longitude location as the number of points a uniform 
grid would have assuming a uniform resolution of AA^. This will produce a uniform grid 
filter function which, through convolution, may be represented by the grid-space weight- 
ing coefficients. The stretched grid weighting stencil may then be constructed by simply 
interpolating the uniform grid stencil to the stretched grid locations. 


The top panel of Fig. 26 shows the grid-space weighting coefficients obtained using this 
empirical method, while the bottom panel shows the differences compared with the eigen- 
value method. The empirical method produces coefficients very similar to the eigenvalue 
method, with the resulting simulation shown in Fig. 27. This implies that the stretched 
grid convolution filter developed in Section 4.1 acts to treat each local A A j as if it were 
global , and constructs a filter function based on a uniform grid of that resolution. Interest- 
ingly the minor differences which are obtained using this empirical method are ultimately 
pathological giving rise to instabilities as shown in Fig. 28 after 30 hours of simulation. 
However, as will be seen in the next section, this technique can be useful to analyze the 
impact of arbitrary uniform grid filter response functions applied to a stretched grid. 


5 Global Shapiro Filter 


As previously noted, the uniform grid GEOS GCM also employs the Shapiro (19 <0) filter 
to globally damp small-scale dispersive waves. This filter also prevents computational non- 
linear instability (Phillips, 1959) to occur. The Shapiro filter is applied as a tendency to 
the winds, potential temperature, and tracers (including specific humidity). Thus, only a 
fraction of the full Shapiro filter is incorporated at each time step. This is done to reduce 
dynamical imbalances and diabatic responses caused by the filter. 
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The Shapiro filter tendency for a quantity q is defined by 


dq 

dt 


SF 



( 36 ) 


Here q and q F are the unfiltered and filtered quantities, and r is an adjustable timescale. For 
the uniform grid GEOS GCM, r is set to 1.5 hours which effectively removes the smallest 
two-grid interval wave in approximately 6 hours. The filter is applied separately in the 
longitudinal and latitudinal direction using 


= [i - i F D n \ [i - (Cx 2 )”] qij . (37) 

where 

{Qij) — ~ — ^qi,j T 9i-l,j) j 

2 1 

(38) 

and n is the Shapiro filter order/2. At 2.5°x2° resolution the GEOS GCM uses an 8th-order 
(n=4) filter. Lower-order filtering corresponds to stronger damping. 


We may examine the Shapiro filter response function in one dimension by assuming a wave 
solution defined by 

q t = q k e lkiAX . (39) 

Doing so we find 

A \ 

= 1 -sin 2 ”(fc— ) . (40) 

The left-hand column of Fig. 29 shows the frequency response for the 8th-order filter on a 
2.5° uniform resolution grid as well as the grid space weights at selected longitudes obtained 
through convolution. We see that the grid space averaging is very local (only 9 grid points) 
and produces a sharp delineation between the small scales which are heavily filtered and 
the longest scales which are virtually untouched. The center column of Fig. 29 shows the 
response of the standard technique described above but applied to the stretched grid. In 
this case the basic local stencil is simply used without concern for the non-uniform mesh. 


Due to the success of the uniform grid interpretation in approximating the convolution 
weights obtained by the eigenvalue method for the stretched grid high-latitude filter, it is 
useful to employ this technique to estimate the proper stretched grid convolution weights 
derived from the uniform grid Shapiro filter frequency response. The right-hand panel of 
Fig. 29 shows the results from assuming that each local AX t was globally uniform. The 
convolution weights associated with the Shapiro filter response for each implied uniform 
resolution were then interpolated to the actual stretched grid locations. The frequency 
response associated with these empirically derived coefficients was then plotted for selected 
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longitude locations. Due to the very local nature of the Shapiro filter and the filter response s 
non-dependency on critical wavelengths, very little impact is seen as a result of the non- 
uniform grid. This is significantly different than the impact seen for the high-latitude filter 
case, Fig. 22, where longitudinal locality greatly influenced filter strength. As a result of this 
analysis in addition to empirical evidence from numerous stretched grid GCM simulations 
and the greater efficiency of the standard algorithm, no obvious advantage is realized by 
using the (empirical) eigenvalue method for the global Shapiro filter. 


6 Conclusion 

In this report the filtering characteristics of a stretched grid GCM have been reviewed. It 
was shown that using the standard FFT polar filtering technique, as performed by Fox- 
Rabinovitz et. al. (1997) for the stretched grid Held-Suarez dynamics benchmark, com- 
putational small-scale noise is generated near the poles in the region of high longitudinal 
resolution. While this noise is associated with high wind speeds in the stratospheric do- 
main, Fox-Rabinovitz et al. found similar problems when total stretching factors were large 
(> 16). For those cases, ad hoc strengthening of the filter coefficients were used to control 
stability. 

A review of the standard high-latitude filtering procedure was made and an examination of 
the physical scales associated with the filter were analyzed on both the uniform grid and 
the stretched grid. It was shown that the implicit assumption of a uniform grid within 
the FFT algorithm results in a misrepresentation of the zonal harmonic composition of the 
input field. Standard filter strengths which are derived from uniform grid considerations 
and stability analyses are thus misapplied to the stretched grid harmonic decomposition, 
resulting in weaker filtering within the high resolution region and stronger filtering in the 
coarse resolution region. This is opposite to the requirement of stronger filtering in areas 
of high resolution. 

A new filtering technique h&s been developed which accurately accounts for the non-uniform 
grid by computing the eigenvectors and eigenfrequencies associated with the stretching. It 
is shown that the convolution of the required damping function within mode space yields 
grid-spaced weights which can be efficiently used to perform the filtering without distortion. 
Offline tests showed that the stretched grid convolution filter correctly filtered more in the 
high resolution region and less in the coarse resolution region. Online GCM experiments 
further showed that computational small-scale noise was no longer generated in the region of 
highest longitudinal resolution near the poles. Away from high latitudes, no adverse effects 
from vising the stretched grid convolution filter were seen. In the limit that the stretched 
grid becomes uniform, this technique reduces to the standard uniform grid filter. 

In addition to the complete eigenvalue/eigenvector solution, an empirical method was de- 
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veloped which treats each local AA* as if it were global The convolution weights associated 
with the filter frequency response for the implied uniform resolution are then interpolated 
to the actual stretched grid locations. This empirical estimation of the convolution weights 
for a stretched grid was used to analyze the impact of applying the Shapiro filter on the 
stretched grid. It was shown that due to the Shapiro filter’s very local grid stencil (9 points 
for the 8th-order filter), the standard application of the Shaprio filter produced very small 
differences compared with the (empirical) eigenvector method. Due to the greater effi- 
ciency of the standard algorithm, no obvious advantage is realized by using the (empirical) 
eigenvector method for global filtering. 
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Appendix 


Linear Algebra Aspects 
of the Stretched Grid Convolution Filter 

The stretched grid convolution filter requires the eigenvalues and eigenvectors of a non- 
symmetric matrix (which we define here as R) from (22). The matrix R is a tridiagonal 
matrix with the addition of the two corner elements rf M and rj~. It has several interesting 
properties. From (19), (20), and (21) it is clear that all row r sums are zero. That is, 
Vz : — 2r; + rf + r~ = 0, or in matrix form, 


Rc = 0 


(41) 


where c is a vector of constant elements c = (c, c, . . - ,c) T . In other words R is singular. 
From the same definitions (19), (20), and (21) it is easily shown the product of the ratios 

rf . ., 

-b is unity: 


n + 

rj 


n 

i= 1 


= 1 


(42) 


since every AA t , AA !+i and AA s _i appears in both the numerator and the denominator. 

If (19), (20), and (21) are all multiplied through by AA;, or equivalently, (22) is multiplied 
on the left by 


| 



" AA X 

1 

A = 

aa 2 

(43) \ 

\ 



i 

i 


a generalized eigenvalue problem containing a symmetric left-hand side R = AR side can 
be obtained: 

R h = A A h . (44) 


It would therefore be straightforward to use a generalized eigen-solver for symmetric ma- 
trices to solve the eigenvalue problem. But the problem can be further simplified. Since 
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both R and A are real and symmetric, it is clear that all A are real. This means that R is 
symmetrizable, that is, a similarity transform can be applied to R to create a symmetric 
matrix with the same eigenvalues: 


Z R Z" 1 = S = S T . 

In order to find Z by construction, consider the diagonal matrix 

Si 

I 

z = 


(45) 


(46) 


Clearly, Z 1 is also a diagonal matrix with j- on the diagonal. 1 If the similarity transform 
is applied using this Z, the result is a matrix with the same structure as R: 


S = Z R Z 



— 2ri 

h r ~ 
S x r 2 

-2r 2 

S 3 V 2 


1 __ 

&1M >r + 

fe r 3 


— 2r IM _j 

Jim r - 


3 

1 



«IM-I IM 


o TM 1 


£im-i 


<>IM 


IM-1 


2r Ti 


(47) 


This matrix is symmetric if, and only if, the elements r+ and r i are non-zero, and the 
following conditions are satisfied: 


«i+i 2 = 



(48) 


The conditions (48) wrap around (i.e., 5\ 2 = ^5 IM ) and can only be satisfied if 

r , 


+ *»+ 


r; r 


1 ' 2 


r, r, 


l ' 2 



(49) 


From (42) we already know this to be the case for the matrix R. One degree of freedom is left 
to set, for example, <*>i = 1. Thus, with a simple similarity transform, a symmetric matrix 

! It is assumed that the are non-zero. 
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can be obtained which is more efficient to factorize 2 than the original R. It also improves 
the algorithm by indicating that the eigenvalues are necessarily real — a non-symmetric 
eigensolver might return very small imaginary components due to round-off errors. Finally, 
the eigen-decomposition of a symmetric matrix yields orthogonal eigenvectors Q, 


Q T S Q = diag(Ai,A 2 ,...,A 1M ) (50) 

which can he a considerable advantage. 

Equation (41) indicates that R has a zero eigenvalue, i.e., Ai = 0. Although most eigen- 
solvers can handle such a case, it is conceivable that floating point errors can result. In 
order to further stabilize the algorithm, it may be worthwhile to remove the zero eigenvalue 
by the following technique. The eigenvector to the zero eigenvalue of R is clearly c, a vector 
of constants, as showm in (41). The equation 

R c = Z” 1 S Zc = 0 (51) 

indicates that y = Zc is the eigenvector of S corresponding to its zero eigenvalue. Assuming 
a normalized y it is possible to create an orthonormal basis, 


Y = [y Y] (52) 

where Y can be constructed iteratively from Gram- Schmidt orthogonalization. Thus a 
second similarity transformation can be applied on S using Y: 


y t sy 


T 


y 

Y 


[s] [y y] 


y T Sy y r SY 
Y r Sy Y r SY 

0 0 T " 

0 s 


(53) 

(54) 

(55) 


Note, for example, that y r SY = y r S r Y = (Sy) r Y = 0 T . This similarity transform en- 
sures that S is also symmetric and contains all the eigenvalues of S except the zero eigenvalue 
corresponding to eigenvector y. The symmetric matrix S has the eigen-decomposition, 


2 For example, with the LAPACK routine dsyev. 


(56) 


Q T S Q = diag(A 2 ,...A rM ) . 

From (53) and (56) it is straightforward to show that the eigenvectors of the non-zero 3 
eigenvalues are YQ. Thus the full set of eigenvectors of S is: 

Q = [ y yq ] . (57) 

It is now possible to specify the eigen-decomposition of R in terms of all the components 
of the fully stabilized eigen-solver: 


M -1 RM 


diag(0, A 

q t sq 

y T 

Q T Y T 


• • • ^im) 

(58) 


(59) 

Z- J RZ[ y YQ ] . 

(60) 


Thus, with the eigenvectors Q of S, the normalized eigenvector y, the set of constructed 

orthonormal vectors Y, and the diagonal matrix Z, it is possible to recreate the eigenvectors 
M of R: 


M = Z [ y YQ J 

which are required in the stretched grid convolution filter. 


(61) 


3 The proof is omitted here that, all other eigenvalues of S are non-zero. 
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Figure 4: Filter function and grid-space weighting coefficients at 85N for uniform grid 
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Uniform— Grid Initial Condition 



Figure 6: Initial conditions for gravity-wave experiment on a uniform grid. 
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Figure 7: Timeseries of single harmonics on a uniform grid using a linearly stable timestep 
and no filtering. 
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Uniform Grid FFT Hlter Lot: 85N 
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Figure 8: Timeseries of single harmonics on a uniform grid using an unstable timestep with 
FFT filtering. 
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Figure 10: A comparison of the uniform and stretched-grid zonal increments. 
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Stretched-Grid Initial Condition 










Figure 12: Timeseries of single harmonics on a stretched grid using a linearly stable timestep 
and no filtering. 
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Stretched Grid FFT Filter Lot: 85N 
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Figure 13: Timeseries of single harmonics on a stretched grid using an unstable timestep 
with FFT filtering. 





Figure 14: Timeseries of single harmonics on a stretched grid using the FFT filter squared. 
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Figure 18: Filter function and grid-space weighting coefficients using standard FFT on a 
stretched grid. 
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Figure 21: Eigenvectors associated with the stretched grid. 
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Stretched Grid Convolution Filter Lot: 85N 
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Stretched Grid Empirical Method Lat 85N 



Figure 27: Timeseries of single harmonics on a stretched grid with weights obtained using 
the empirical method. 
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Eigenvalue Method Wavenumber 1 


Empirical Method Wavenumber 1 



Eigenvalue Method Wavenumber 3 



Eigenvalue Method Wavenumber 5 



Eigenvalue Method Wavenumber 7 
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Figure 28: Single harmonics at time = 30 hrs comparing the eigenvalue and empirial meth 












Figure 29: Filter function and grid space weights for an 8th-order Shapiro filter on the 
uniform and stretched grid. 
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